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For  a  depth  dependent  reference  speed,  c(z),  we  can  no  longer  invert 
the  integral  equation  exactly.  However,  we  can  write  down  an  asyaptotic 
i  high  frequency  approxiaation  for  the  kernel  of  the  integral  equation  and  an 

asyaptotic  solution  for  the  perturbation.  The  coaputer  iapleaentat ion  of 
this  result  is  designed  along  the  saae  lines  as  the  code  for  the  constant 
background  case.  In  tests  of  processing  tine,  we  find  that,  at  worst,  the 
total  processing  tiae  for  this  algoritha  with  depth  dependent  background 
soundspeed  is  about  the  saae  as  for  a  coaparably  prograaaed  k-f  algoritha 
with  constant  background.  By  worst  we  nean  that  we  choose  the  aspect  ratio 
—  the  nuaber  of  traces  divided  by  the  nuaber  of  points  per  trace  —  to  be 
optiaal  for  the  k-f  algoritha.  Ve  present  exaaples  which  deaonstrate  the 
aethod  iapleaented  as  a  aigration  technique  and  coapare  with  the  application 
of  alternative  aigration  algorithas.  The  exaaples  we  chose  were  ones  in 
which  the  objective  is  to  iaage  the  flanks  of  a  salt  doae. 


— — /'The  purpose  of  this  paper  is  to  describe  an  ezteasioa  of  the 
aultidiaensionl  Born  inversion  technique  [Cohen  and  Bleisteia,  1979a]  for 
acoustic  waves.  In  that  earlier  work,  a  perturbation  in  reference 
soundspeed  was  deterained  by  scanning  that  the  reference  or  background  apeed 
was  constant.  In  this  extension,  we  allow  the  reference  apeed  to  be  a 
function  of  the  depth  variable,  x,  bnt  atill  require  that  it  be  independent 
of  the  transverse  variables.  The  ontpnt  of  this  aethod  is  a  high  frequency 
bandliaited  reflectivity  function  of  the  subsurface.  The  reflectivity 
function  is  an  array  of  bandliaited  singular  functions  sealed  by  the  noraal 
reflection  strength.  Each  singular  function  is  a  Dirac  delta  function  of  a 
scalar  arguaent  which  aeasurea  distance  normal  to  a  reflecting  interface. 
Thus,  the  reflectivity  function  is  an  indicator  asp  of  subsurface  reflectors 
which  is  equivalent  to  the  asp  produoed  by  migration.  In  addition  to  the 
assumption  of  aaall  perturbation,  the  aethod  requries  the  aasuaption  that 
the  reflection  data  reside  in  the  high  frequency  regiae,  in  a  well-defined 
tenu^ 

The  aethod  is  based  on  the  derivation  of  an  integral  equation  for  the 
perturbation  in  aouadapeed  froa  a  known  reference  apeed.  Then  the  reference 
speed  is  constant ,  the  integral  equation  adaits  an  analytic  solution  us  a 
aultifold  integral  of  the  reflection  data.  Further  high  freqneney 
aayaptotie  analysis  simplifies  this  integral  considerably  and  leads  to  an 
extreaely  efficient  nuaerieal  algoritha  for  eoapnting  the  relfeetivity 
function.  In  a  paper  by  Blaisteia,  Cohen  and  login  (19141.  the  development 
of  a  computer  code  to  iapleaent  this  constant  reference  speed  solution  is 
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1.  IN1MMCTI0N 


The  purpose  of  this  psper  is  to  describe  sn  extension  of  the  Born 
inversion  technique  [Cohen  end  Bleistein,  1979a].  In  that  paper,  a  nethod 
was  developed  to  deteraine  the  perturbation  in  soundspeed  f  ron  a  known 
reference  speed.  A  closed  fora  solution  was  derived  when  the  reference 
speed  was  taken  to  be  constant.  In  this  extension,  we  allow  the  reference 
speed  to  be  a  function  of  the  depth  variable  z,  but  we  require  that  it  be 
independent  of  the  transverse  variables  (laterally  hoaogeneous) .  The 
perturbation  is  allowed  to  be  a  function  of  all  three  spatial  variables  when 
a  planar  survey  on  the  upper  surface  is  available i  it  is  allowed  to  be  a 
function  of  one  transverse  variable  and  depth  when  only  a  line  of  surface 
data  is  available.  Ve  shall  refer  to  this  latter  case  of  three  diaensional 
propagation  over  an  earth  with  two  diaensional  variations  as  the  two  and- 
one-half  diaensional  eases  we  shall  refer  to  the  foraer  as  the  three 
diaensional  case.  In  both  cases,  the  surface  observations  we  use  are 
backseattered  (CMP  stacked)  tiae  logs. 

Let  us  denote  the  sonadspeed  by  v(x,y,z).  We  introduce  the  reference 
speed,  c,  and  a  perturbation,  a,  in  the  fora 

-L  -  -L  (1  +  «)  (1) 

v  e 

The  aain  features  of  the  aethod  are  as  follows: 

(i)  The  aethod  starts  froa  the  wave  equation  and  a  nonlinear  integral 
equation  for  o(x,y,s)  is  derived.  The  equation  is  nonlinear 


because  it  contains  the  product  of  the  perturbation  and  the 
unknown  interior  wave  field. 

(ii)  Linearization,  under  the  assumption  that  a  is  "small,*  leads  to  an 
integral  equation  for  a  in  which  the  data  are  the  ensemble  of 
backscatter  (coincident  source-receiver)  observations  on  the  upper 
surface.  The  kernel  of  the  integral  equation  is  the  square  of  a 
Green's  function  for  the  unperturbed  problem  with  reference  speed 
c(z) . 

(iii)  For  constant  c.  the  integral  equation  admits  a  closed  form 
solution  for  o  as  a  multi-fold  intefral  of  the  observed  data 
[Cohen  and  Bleistein,  1979].  Asymptotic  analysis  [Bleistein, 
Cohen  and  Hagin,  1984*  referred  to  as  BCH,  below]  reduces  the 
number  of  integrations  which  must  be  performed  by  computer  so  that 
the  processing  times  are  competitive  with,  or  better  than,  times 
for  existing  migration  codes.  (See  Gray,  [1984].)  For  depth- 
dependent  c,  the  integral  equation  admits  an  approximate 
(asymptotic)  solution,  again  as  a  multi-fold  integral  of  the  data. 
As  with  the  case  of  constant  c,  asymptotic  analysis  greatly 
simplifies  the  computer  implementation. 

(iv)  The  computer  i^lementation  takes  account  of  the  fact  that  real 
seismic  data  is  bandlimited  in  frequency.  The  method  also  employs 
constraints  which  aeeouat  for  causality  and  avoid  aliasing  in  the 


The  output  of  this  method  is  a  high  frequency  bandlimited  reflectivity 


fraction  of  the  nkmfis*.  Th#  wfUetitlty  fraction  is  an  array  of 
baadliaited  sinanlar  fractions  scaled  by  the  noraal  incidence  reflection 
strength.  Bach  singular  fraction  ia  a  Dirac  delta  fraction  of  a  scalar 
argraent  which  aeaenres  distance  noraal  to  a  reflecting  interfaces  see  Fig. 
1. 


In  characterising  this  data  aa  high  freqneney  data  we  naan  that  the 
diaeneionleaa  paraaeter  ■- 


X  ■  4xfL/c 


should  be  large.  In  this  equation,  f  ia  the  ainiara  frequency  in  Is,  o  it 
the  local  aoradepeed  and  L  is  a  typical  length  scale.  There  is  an  extra 
factor  of  2  in  this  equation  due  to  the  two-way  travel  tiae  of  interest  in 
inverse  pr obi east  alternatively,  c/2  nay  be  viewed  as  the  "nigra t ion 
velocity"  for  a  problea  where  sources  are  set  off  at  each  reflector  at  tiae 
sero.  The  length  sealea  of  interest  are  the  range  froa  the  upper  surface  to 
the  reflecting  surface,  the  principal  radii  of  curvature  of  the  reflector  or 
the  distance  between  reflectors.  For  the  first  two  of  these,  a  lower  liait 
of  SOOsi  is  reasonable.  Let  us  consider  a  lover  liait  of  frequeaoy  of  6  Is 
and  a  aoradepeed  of  2000  a/ see.  For  this  choice,  X  ■  6s  is  aneh  larger  than 
necessary  for  high  frequency  ssyaptotie  aethods  to  be  valid.  In  feet,  a 
value  of  three  (or  s)  will  usually  suffiee.  Thus,  even  an  increase  of 
velocity  by  a  factor  of  S  or  a  decrease  in  length  scale  by  the  suae  aaorat 
would  leave  X  large  eaeugh  for  ssyaptotie  aethods  to  apply.  Unfortunately, 


11  offset  seisaie  data  do  not  usually  support  resolutiea  of  layers  whose 


t 


separation  is  such  as  to  sake  X  less  than  three,  whether  or  not  we  exploit 

high  frequency  techniqnes  to  invert  the  data.  - 

i 

The  aethod  has  been  checked  on  synthetic  data  in  the  two  and-one-half 
dimensional  case.  The  ontpnt  of  this  inversion  algorithm  is  a  depth  profile 
equivalent  to  one  produced  by  a  depth  migration.  The  added  feature  of 

inversion  is  that  the  amplitude  of  the  output  provides  an  estimate  of  the 
velocity  increments  across  the  interfaces.  The  accuracy  of  this  estimate 
can  be  no  better  than  the  accuracy  of  the  input  data  as  relative  true 

amplitude  data.  However,  in  BCH,  we  show  theoretically,  using  asymptotic 
analysis,  that  for  true  amplitude  data  the  output  provides  a  more  accurate 
estimate  of  reflection  strength  at  interfaces  than  its  basis  in  perturbation 
theory  would  suggest.  In  the  absence  of  true  amplitudes,  the  output  is 
fully  equivalent  to  the  output  of  a  migration  algorithm.  In  BCH,  we  discuss 
the  relationship  between  inversion  of  both  k-f  migration  [Stolt,  1978]  and 
Kirchhoff  migration  [Schneider,  1978].  In  order  to  distinguish  the  case  in 
which  parameter  estimation  is  possible  from  the  case  in  which  it  is  not,  we 
shall  refer  to  the  former  as  a  seismic  inversion  and  the  latter  as  a 

structural  inversion.  Thus,  from  our  point  of  view,  migration  is  a 
structural  inversion  technique. 

For  a  depth  dependent  reference  speed,  c(z),  we  cannot  invert  the 
linearized  integral  equation  exactly  for  u.  Ve  cannot  even  write  down  an 
exact  analytical  expression  for  the  kernel  of  the  integral  equation,  which 
is  proportional  to  the  square  of  a  Green's  function  for  the  wave  operator 

i 

with  sowadspeed  c(z).  However,  we  can  write  down  an  asymptotic  (high 


inversion  of  the  integral  equation  for  a(x,y,z).  Indeed,  fnndaaental  to 
this  extension  is  early  exploitation  of  asyaptotics  on  the  basis  of  high 
frequency  in  fall  anticipation  of  obtaining  a  solution  only  in  the  high 

frequency  regiae  and  only  for  the  bandliaited  reflectivity  function.  (This 
is  auch  in  the  spirit  of  Clayton  and  Stolt  [1981]).  When  the  data  is  not 
relative  true  aaplitode  data,  we  lose  only  the  estiaate  of  reflection 
strength  and  this  seisaic  inversion  aethod  reaains  a  structural  inversion 
aethod,  now  in  the  context  of  a  depth  dependent  background  soundapeed. 

The  coaputer  implementation  of  this  result  is  designed  along  the  saae 
general  lines  as  the  constant  reference  speed  algorithm  described  in  BCH. 
To  obtain  the  output  at  a  subsurf aoe  point,  it  is  necessary  to  take  the 
temporal  Fourier  transform  of  the  range  normalised  data,  multiply  by  a 
filter  deduced  by  the  theory,  invert  the  transform,  evaluate  each  trace  at  a 
travel  time  deduced  from  the  theory  with  an  appropriate  amplitude  scale,  and 
integrate  over  a  set  of  traces.  In  the  constant  reference  speed  case,  the 
traveltime  and  amplitude  scales  of  the  last  integration  are  explicit 
functions  of  the  integration  point  and  the  outpnt  point.  In  the  e(s) 

reference  speed  case,  these  latter  functions  are  replaced  by  the  geometrical 
optica  traveltime  and  an  amplitude  on  the  ray  connecting  the  subsurfaoe 

outpnt  point  and  the  aonroe/reeeiver  point  in  the  context  of  the  depth 

dependent  reference  profile.  These  functions  are  given  implicitly  in  terms 
of  a  (ray)  parameter  which  remains  constant  on  a  ray  and  is  determined  from 
the  equation  of  the  ray  connecting  the  two  points.  Ve  arbitrarily 
characterise  the  reference  speed  in  this  iavleaentation  as  being  pieeeviee 
linear  and  continuous,  connecting  prescribed  values  of  souadspeed  nt 
prescribed  depths  in  the  subsurfsce.  (For  structural  inversion,  we  could 
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just  at  easily  prescribe  the  reference  velocity  to  be  piecewise  constant. 
For  seisaic  inversion,  the  piecewise  constant  background  would  required  a 
mote  coaplicated  aaplitude  in  the  inversion  operator.)  Ve  coapute  the  ray 
paraaeter  and  the  necessary  functions  in  advance  and  retain  then  in  tables 
to  be  called  as  needed  for  the  final  integration  of  the  algoritha. 

Because  these  tables  need  only  be  coaputed  once  for  each  output  depth 
(i.e.,  for  each  transverse  offset  separating  output  point  and 
source/receiver  location),  the  coaputing  tiae  for  this  part  of  the 
processing  becoaes  a  progressively  saaller  fraction  of  the  total  processing 
tiae  as  the  aspect  ratio  —  the  nuaber  of  traces  divided  by  the  nuaber  of 
output  points  —  becoaes  larger.  In  tests  of  processing  tiae  we  find  that, 
at  worat,  the  total  processing  tiae  for  this  algoritha  with  depth  dependent 
background  velocity  is  about  the  saae  as  for  a  coaparably  prograaaed  k-f 
algoritha  with  constant  background  velocity.  By  worst,  we  Bean  that  we 
choose  both  the  nuaber  of  traces  and  the  nuaber  of  points  per  trace  to  be  a 
power  of  two  and  we  choose  the  aspect  ratio  —  nuaber  of  traces  divided  by 
nuaber  of  points  per  trace  —  to  be  optiaal  for  the  k-f  algoritha. 

Because  of  its  iapleaentation,  described  above,  the  awthod  bears  close 
relationship  to  one  described  by  Carter  and  Frazer  11984].  The  chief 
difference  between  then,  other  than  algorithaic  details,  are:  (1)  Carter 
and  Frazer's  aethods  allows  for  transverse  variations  in  c  while  ours  does 
nott  (2)  our  aethod  allows  for  processing  of  true-aaplitude  data  for 
reflection  strength  and  henoe  velocity  while  Carter  and  Frazer's  does  not. 

We  present  synthetic  data  ezaaples  which  deaonstrate  the  aethod 
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implemented  as  a  structural  inversion  (migration)  technique.  For  the 
simplest  example  of  a  single  steeply  tilted  reflector,  ve  compare  this 
method  with  output  of  a  k-f  algorithm  applied  to  a  time-stretched  version  of 
the  same  data  and  to  a  Eirchhoff  migration.  Ve  alao  present  two  examples  in 
which  the  objective  is  to  image  the  flanks  of  a  salt  dome  intruding  into  an 
otherwise  horixontally  stratified  medium.  The  advantagea  of  the  method  of 
this  type  of  application  are  evident  from  the  output. 


2.  MUTATION  OF  TIE  IHIKKAL  EQUATION  FOB  ALPHA 


We  shall  derive  aa  integral  equation  for  a.  We  atsaae  that  the  tiae- 
rednced  wavefield  u(x.  »)  if  a  solution  of  the  Belaboltz  equation 


* 

V  *u  +  ~  n  *  -6(x  -  $)  6(y  -  t»)  A(x) 
v 


(2) 


In  this  equation,  V*  is  the  Laplace  operator  and  6  is  the  Dirac  delta 
function.  This  equation  is  to  hold  for  all  (x,y,z).  It  is  further  aasuaed 
that 
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v*  c*(0) 


x  >  0,  a  **  a(x,y.z),  c  “  c(x)  i 

x  <  0  . 


(3) 


The  function  c(x)  is  the  known  continuous  reference  velocity  and 
«(x,y,x)  is  the  unknown  perturbation  to  be  deterained  from  observations  on 
the  upper  surface. 

We  introduce  uj»  the  solution  of  the  unperturbed  Relaholts  equation 
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The  function  Uj  i«  often  referred  to  ee  the  "incident  me.*  Then  we 
■ay  write 


a  -  Oj  +  ag  » 


(5a) 
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+  a. 
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Let  aa  sappose  that  a  ia  aaali.  Hear!  at  lea  lly,  aiaee  the 
"aooree  ten”  -(*»*/c*)oa  ia  eqoatioa  (5b)  ia  0(a),  we  expeot  that  the 
aolatioa  eg  ia  aa  well.  Thai  aa  a  firat  approxiaation,  we  aay  aegleet  the 
prodact  oag  oa  the  right  aide  aa  being  qaadratio  ia  ■  aad  heaee  of  lower 
order  than  the  prodaet  oaj.  Conaeqaaatly,  the  liaearixed  eqoatioa  for  «g  i* 


(«) 


Ve  ahall  bow  write  dowa  the  Oreea'a  faaetioa  repreaeatatioa  for  the 
aolatioa  ig  eralaated  book  at  the  aoaree  point i  thia  ia  the  baekaeatter 
reapoaae  which  ia  aaaaaed  to  be  the  obaerwed  data.  Ve  reaark  that,  frea 
(2),  the  Oreea'a  faaetioa  ia  jaat  aj  iteelf.  aad  the  repreaeatatioa  for  ag 


it 
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“  -***  |  dx  |  dy  |  di  gj(i,y,zi(,q>Oit) 


(7) 


In  this  iqwtlra,  w  have  nt«4  the  feet  that  a  it  nonxeto  for  t  <  0  to 
integrate  only  over  positive  z  values.  Alto,  we  have  exhibited  both  sett  of 
spatial  argnaents  of  ij  to  eaphasixe  its  dependence  npon  both  observation 
and  integration  coordinates.  This  is  the  integral  equation  for  ai  it  is  the 
first  tern  in  the  Born  series  for  Ug.  We  propose  to  solve  this  integral 
equation  in  the  high  frequency  regiae. 


10 


3.  ANALYSIS  OP  UK  INIB6KAL  EQUATION  PON  ALMA 


Aa  exact  analytical  inveraion  of  Eq.  (7)  for  o  ia  not  generally 
available  unleae  c  ia  conataat.  Thua.  we  aoat  ooatent  onraelvea  with  aa 
approxiaate  eolation  which  retaine  featnrea  of  intereat  of  the  true  veloeity 
profile.  Aa  diacnaaed  in  the  introduction,  we  ahall  develop  a  high 
frequency  aolutlon  for  the  reflectivity  function  under  the  aaauaption  of  a 
piecewiae  conatant  a. 

Ve  envieion  an  inveraion  algoritha  whioh.  for  output  a(x.y.x).  will  be 
in  the  fora  of  an  integration  over  a  aet  of  traeea  whoae  traaaverae  liaita 
define  aoae  bounded  interval  about  the  traaaverae  coordinate  of  the  output 
or  obaervation  point  -  a  Kirchhoff  integral.  Toward  thia  end.  we  ahall  be 
content  to  replace  uj  in  (7)  with  a  high  frequency  baadliaited  approxiaatioa 
of  that  function  for  which  the  tranaverae  offaet  ia  bounded.  Thia 
approxiaatioa  ia  the  downward  propagating  WKB  or  geoaetrioal  optica 
approxiaate  aolutlon. 

In  Appendix  A  we  derive  aueh  a  aolutlon  aa  a  Fourier  integral.  He 
reault  ia 


uT(x.y.ii{.i|»Oi«)  ■  ~  — ~ 
1  4**i 


dp.  dp,  exp  {ip* (x  -  l)  ♦  i#(x,pt«)) 

— " —  ; .  ■  ■  ■_  - -  •  <•> 

^P,<*#P»*>  "p,W»*»*> 


-II- 


Ill  this  equation,  we  have  introduced  the  following  notation: 


S  -  U 


agn  w 


^0 


da'  . 


|i((z,g,w)  - 


w _  * 

1 - H  * 

c  (a) 


i  >  i 
U  “  +  * 


(9) 


Ve  do  not  concern  ourselves  hare  with  the  values  of  jt  for  whieh  the 
square  root  is  zero  or  imaginary.  The  former  ease  corresponds  to  rays  at 
their  turning  point t  the  latter  ease  corresponds  to  evanescent  waves.  Ve 
justify  the  first  of  these  by  the  observation  that,  in  our  inversion 
formula,  we  shall  artificially  limit  the  range  of  transverse  variables  to 
preelude  the  turned  rays,  thus  limiting  the  amount  of  dip  whieh  can  be 
imaged  to  vertieal,  at  most.  Ve  justify  the  latter  by  noting  that,  in  fact, 
«  is  the  Fourier  inversion  of  its  spatial  transform  over  real  wave  vectors. 


By  using  this  representation  for  uj  iu  (9),  we  obtain  the  following 
approximate  integral  equation  for  u: 


-U- 


dx  djr 


dx 


Bg(£.0|«*) 


4(2x) 


I  c*(x) 
J0 


d|it  d(i4  exp{ig’(x  -  £)  +  id(x>||t«)} 


(10) 


dXx  dX,  exp{i£*(x  -  i>  +  i#(x, £»•»)) 

R,(0,&»o> 


la  till  equation,  ve  have  omitted  Units  of  integration  exoept  for  the  loner 
limit  in  x.  We  take  the  point  of  view  that  the  limits  of  integration  in  all 
integrals  is  the  range  of  valwi  that  allows  the  square  roots  in  the 
integrand  to  remain  real. 

Eq.  (10)  is  the  integral  equation  which  we  shall  invert  in  the  next 
section.  Ve  note  that  Bq.  (10)  is  a  'high  frequency*  approximation  to  Bq* 
(7),  which  is  itself  a  "small  o*  approximation  to  the  exact  expression  for 
«g((,i|,0ie).  Ve  also  note  that  we  shall  invert  in  Bq.  (10)  only  in  an 
asymptotic  sense.  The  validity  of  our  cascade  of  approximations  will  ho 
demonstrated  hy  the  numerical  examples. 


DRBRAL  BQDATION 


We  shall  now  describe  the  (asyaptotic)  invert  ion  of  the  integral 
equation  (10).  To  begin,  we  introduee  the  apatial  Fourier  transform  of  ag 
with  reapeet  to  J  defined  by 


Vg(k.O)M) 


|  d(  dq  Ug(^.OjM)  exp[-2ik’£)  ,  k  *  (k^.k^) 


(11) 


Bw  taking  this  Fourier  tranafom  in  (10)  we  obtain  the  following  equation: 


«s(k.0i«t) 


m  m  m 

*  [  I  •  f  a(x,z) 

-r-  I  **  I  M*  — r — 

*  C  (*) 
J0 

1 


dpt  dp,  exp(i||*x  4  i#(*,giw)) 
^|i,(x,lti«)  p|(0,|ttw) 


1dXx  dX#  ezpU&’g  4  ip(x,^»t»)) 
^l»» <**&»ee>  M,(0,£iw> 


(12) 


M2fc  +  h  +  tt>  • 


Iu  this  result,  we  have  used  the  feet  that 


The  delta  function  ia  (12)  now  allows  as  to  oalcalate  the  integral  ia  V. 
(We  eoald  as  eaailjr  aae  it  to  oalealata  the  iategral  ia  g>)  The  resnlt  is 


We  shall  oalealata  the  g  iategral  by  the  method  of  stationary  phase 
[Bleiateia,  1984].  The  large  parameter  in  this  iategral  ia  agaia  the  high 
frequency.  although  that  ia  aot  apparent  ia  the  given  representation.  It 
would  beeoae  apparent  if  we  were  to  seals  g  and  |  by  I  eel  and  seals  the  depth 
by  the  depth  to  the  first  refleetor.  Then  we  would  find  the  dimensionless 
depth  to  the  first  reflector,  measured  ia  wave  lengths*  emerging  as  a  scale 
oa  the  phase  of  the  g  integration  and  this  would  be  the  dimensionless  large 
parameter.  We  do  aot  do  this*  but  proceed  formally  to  use  the  method  of 
stationary  phase  directly  on  the  iategral  (1$).  This  ia  carried  out  in 
Appendix  B.  The  result  ia 


Og(k » 0 ; w) 


dz  8(k,z)  exp  (  2ifKz,kt«»)  } 


(14) 


o> 

ftni 


c  (z)  |i((z.k;u)  ^(O.kjw)  o(z,ktw) 


In  thiz  equation. 


(IS) 


When  c(z)  it  replaced  by  a  constant,  the  exponent  in  tbe  integrand  in 
(14)  is  siaply  ik,z,  with 


k 


s 

i 


Tbe  integrand  also  siaplifiea  significantly  is  this  ease.  One  then 
recognises  (14)  as  an  eqnation  for  the  Ponrier  t  tans  fora  with  respect  to  a 
of  V.  The  inversion  of  this  transfora  to  yield  n  is  then  straightforward. 
Awareness  of  this  fast  snggasts  an  inversion  which,  st  least  asymptotically, 
has  the  same  affect  when  a  varies  with  depth.  This,  we  define 


w(£*t> 


(16) 


16i  d*k  d"  exp{-2i#(f,|jt*)  +  2ik*5) 

n*  4,(t.k»“>  sgn  * 

We  look  upon  v  at  the  result  of  applying  an  integral  operator  to  the 
(Fonrier  transformed)  data  in  (14).  By  applying  the  same  operator  to  the 
right  side  of  (14).  we  find  that 


•  •  • 

"  V  I d*  I dy  I  d“  I dkl  I dk*  I  «*<*>  n,<t.s 

*—m  -•  *0  J  J  J 


«(*,*)  exp(2i[#(s.ki»)  -  p(f.kiw)]  +  2ik*({  -  j)} 
j»1(z.kj«)  |i ,(0,kj»»)  o(a,|»»> 


For  the  eaae  of  constant  e.  the  integrations  on  the  right  aide  ean  he 
carried  out  explicitly  to  yield  the  result 

w(J.f)  ■  c  a(£,f)/ft  a(£,f)  ■  t  w(£»t)^«  •  (*•) 

Thus  at  least  for  constant  c.  we  see  that  the  multifold  integration  on  the 
data  represented  by  w  is  in  fact  a  multiple  of  at  that  is,  (16)  and  (IS) 
prowide  a  solution  for  a. 

Our  object iwc  now  is  to  carry  out  the  integrations  in  (17)  for  the  ease 
•f  variable  e.  Again,  we  do  so  by  using  the  method  of  stationary  phase. 
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this  tiae  in  the  four  vsrisbles  k  and  The  reaaining  integrals  in  <0  and  z 
then  becoae  straightforward.  This  analysis  is  carried  out  in  Appendix  C. 
The  result  is 


w<j.r> 


C (f)  c(0)  0(4, t> 


c(z)  dz 


(19) 


which  provides  the  solution  for  a: 


m{i’V  m  cTT)  Ci6)  c(,)  d* 


(20) 


In  order  to  write  the  solution  for  a  in  tens  of  the  data  observed  at 
the  upper  surface,  we  use  this  result  and  the  definition  of  w  in  (16)  along 
with  the  definition  of  V  in  (11)  and  the  definition  of  the  tiae  transfon 
which  placed  us  in  frequency  doaain  to  begin  with.  The  result  is 


l<i  4  •<*>  f  f  f  C  f  C 
c(f>  c(0)  J  1  J  *  J  L  L  J0 


dt 


(21) 


BjCj.O.t)  exp(2itk’(£  -  -  aga  a  Jq  y, (*,&•«> 


•*»  •  !*,<£»!»•> 
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Ia  this  equation,  we  have  ssed  the  notation  Og  to  denote  the  upward 
scattered  field  in  the  tine  domain. 


As  aentioned  in  the  introduction,  we  do  not  process  the  data  to  yield  a 
itself,  but  to  produce  the  reflectivity  function  or  normal  derivative  of  a, 
at  each  of  the  interfaces  between  the  layers  of  constant  a.  The  theory 
developed  in  Cohen  and  Bleisteia  11979b]  shows  that  we  can  processs  for 
da/d n  by  aultiplying  the  integrand  by  -21s/«(f)i  we  can  obtain  3  if, 
instead, we  multiply  by  -im/2e(f)  before  inverting.  This  is  also  discussed 
in  BCH  and  in  Bleisteia  [1984].  Thus,  we  introduce  the  reflectivity 
function  obtained  by  aultiplying  the  integrand  in  (21)  by  this 

factor.  That  is. 


Mtr>  ■  L‘y  I? 

ft 

V*.O.t>  e*p(2i[k* (£  -  *)  -  sga  a  JQ  +  i**t) 


(22) 


*8®  • 


In  practice,  an  areal  array  of  data,  as  is  required  for  this  formula, 
is  not  usually  available  i  only  a  line  of  data  at  the  upper  surf  see  is 
available,  aay  ia  the  direetioa  of  the  s  axis  of  our  formulation.  In  order 
to  treat  this  type  of  data,  it  ia  assumed  that  there  is  no  variation  in 
velocity  ia  the  orthogonal  horisoatal  direetioa,  namely  the  y  direetioa.  In 
that  ease,  the  data  itself  asst  be  y- independent.  This  allows  us  to  carry 
out  the  integrations  in  y  and  k,.  The  former  of  these  yields  *  b(h,>  »•* 
the  latter  integration  then  yields  a  value  of  n  while  also  evaluating  the 
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iat«|md  it  k,  -  0.  Consequently,  «•  ujr  npliei  k1  by  k  in  the  integrsnd 
•ad  write  dowa  the  foilowiag  result  for  the  refleetivity  faaetioa  ia  the 
eeee  of  y  iadepeadeace : 


S  o(z)  ds 

a  e*(t>  d(0) 


Ug(x,0,t)  exp(2i[k((  -  x)  -  sgnw  JQ|if  (z,ki»»)  dxl  +  iat) 


sgn  w  |*t  (z.ktw) 


(23) 


Although  9  is  sssuaed  only  to  have  two  diaeasioael  variability  ia  this 
result,  the  wave  field  has  three  diaeasioaal  spreadiag  properties.  Thus, 
this  is  the  foraula  for  the  two  aad  one-half  diaeasioaal  case.  Te  reaark 
that  a  fully  two  diaeasioaal  foraalisa  (goveraed  by  a  two  diaeasioaal  wave 
equation)  would  not  properly  seeouat  for  the  "geoaetrical  spreadiag*  of  a 
three  diaeasioaal  eaviroaaeat  as  this  two  aad  oae-half  diaeasioaal  result 
does. 


Por  the  ease  of  coastaat  o,  we  esa  coap are  this  foraula  to  the  result 
ia  Cohea  aad  Bleisteia  [1979a].  This  eoaparisoa  is  suffieieatly  eoaplieated 
to  carry  out  that  we  do  so  ia  Appendix  0.  To  leading  order  asyaptotioally, 
they  do,  indeed,  agree. 

Also  for  constant  o,  (23)  has  the  fora  of  a  k-f  aigration  integral.  Ia 
fact,  when  eoapared  to  ftolt  [1978],  equation  (SO),  we  find  that  we  auat 


(i)  interpret  Stolt't  ?  at  tUg(x,0,t)  end  <ii)  integrate  by  perte  with 
respect  to  a  ea  it  done  in  Appendix  D  of  this  paper.  Then  the  integrals  are 
the  sane  to  leading  order  asymptotically  for  large  w  up  to  a  multiplicative 
constant.  Since  the  Stolt  k-f  fonsnla  it  aigrating  the  field  and  wa  are 
producing  0,  this  difference  of  a  multiplier  censes  no  difficulty,  however, 
except  for  such  constant  aultipliers,  it  is  necessary  that  the  integral 
processing  agree,  at  least  to  leading  order  in  s,  because  both  methods 
produce  an  image  of  the  interfaces  between  layers.  Only  the  interpretation 
of  the  intensity  of  the  iaaging  function  differs  in  the  two  aethoda. 
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The  font  via  (29)  is  not  easily  implemented  on  the  eompnter. 
Putkeraore,  ia  that  formula,  we  have  aot  tally  exploited  the  poeaible 
aimplif ieatioaa  available  aader  the  assumption  of  "high  frequency."  To 
aehieve  this  end,  there  ia  still  one  integration  which  can  be  earried  out 
asymptotically,  namely,  the  integration  in  k,  the  transverse  wave  number. 
This  calculation  is  carried  out  ia  Appendix  B.  The  result  of  that 
calculation  is 


32  Jgc(s)  dx 
f*  c(0)  c‘<t) 


dx 

£*<t)  -  p*  8<t.P> 


|jR 

|  Ug(x.O,t) 


dm  exp(~2imr(p,f)  +  i(sgnm)x/4) 


exp(imt)dt  i 


(24) 


I*  -  Cl 


8<C.P> 


£ 


(»)  dx 

•  .3/2 


(n  (a)  -  p  ) 


1/2 


.»(,)  -£<2>. 

C  (l) 


(25) 


In  these  equations,  the  integrand  is  stated  parametrically  with  parasMter 
p.  This  is  a  dimensionless  ray  parameter  which  denotes  the  sine  of  the 
emergence  angle  with  respect  to  the  vertical  of  a  reflected  ray.  (This 
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usage  is  a  departure  from  the  usage  in  the  t-p  method  where  this  parameter 
is  divided  by  c(0)  but  is  also  denoted  by  p.)  For  each  choice  of  |x  -  (| 
■ad  we  mutt  determine  the  value  of  p  from  the  first  equation  in  (2S). 
Given  p,  the  values  of  t(t,p)  and  S(f,p)  are  then  determined,  as  well  as  the 
explicit  square  root  appearing  in  the  integrand.  The  choice  of  t  determines 
a  "time"  at  which  the  processed  time  logs  (CDP  stacked  traces)  are  to  be 
evaluated.  'Processed*  means  Fourier  transformed  in  time,  multiplied  by  a 
function  of  o>  and  inverse  Fourier  transformed.  The  range  of  integration  in 
space  is  restricted  to  those  values  of  x  for  which  the  square  roots  of  the 
integrand  remain  real  and  nonxero  and  the  travel  time  z  remaining  leas  than 
the  maximum  time  on  the  set  of  traces. 

Equations  (24)  and  (25)  are  to  be  implemented  on  the  counter.  In 
anticipation  of  this,  we  have  written  our  result  in  terms  of  n(x),  the  index 
of  refraction,  which  is  a  number  varying  typically  between  .25  and  4.  rather 
than  in  terms  of  the  inverse  velocity  or  slowness,  which  is  three  orders  of 
magnitude  smaller,  with  similar  estimates  true  for  the  corresponding 
parameter  p.  The  computer  imp 1 eme ntation  proceeds  as  follows: 

(i)  Calculate  the  Fourier  time  transform  of  each  traoe. 

(ii)  Multiply  by  a  function  of  the  frequency  variable,  perform 
filtering  as  required  and  take  the  inverae  Fourier  tranaform. 
This  inverae  transform  is  now  given  on  a  uniform  grid  in  time. 

(iii)  For  fixed  output  depth  f,  develop  a  table  of  values  of  p  as  a 
function  of  |s  -  (|  (transvexae  offset  between  output  point  ({.f) 
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•ad  neiifir  loeation  (x,0))  for  the  prescribed  rtftrnct  velocity 
c(x).  Ve  ace  a  table  here  becaaae  the  saae  choice*  of  there 
variable*  arise  repeatedly  in  carrying  oat  the  spatial  integration 
for  different  oatpat  points  at  the  given  depths. 

(iv)  With  the  p  valaes  deterained,  develop  tables  of  the  other 
faaetions  of  p,  namely  t  and  S  in  Eq.  (25),  aa  well. 

(v)  Nov  for  fixed  ($,f),  calculate  the  spatial  integral  a*  a  discrete 
san  over  the  locations  x  of  the  traces,  sabject  to  the  liaits  of 
saaaMtion  noted  above.  Typically,  the  valae  of  x  at  a  given 
choice  of  x,  (  and  f  will  not  be  available  in  the  table  of 
processed  data  deterained  in  (ii).  Therefore,  nse  an 

interpolation  scheae  to  estimate  the  valae  of  the  processed  data 
at  i  in  terms  of  the  processed  data  at  the  grid  points  of  the  time 
variable. 

The  integration  in  e  is  to  be  carried  oat  over  the  available  baadvidth 
of  the  data.  The  Foarier  transform  shoal d  not  simply  be  traacated  at  the 
endpointa  of  the  bandwidth,  bat  a  smoothing  filter  shoald  be  applied  to 
•void  ringing  of  the  oatpat  dae  to  a  discontinaoas  filter. 

The  spatial  integration  aaat  also  be  farther  constrained  by  sampling 
rate  considerations  associated  with  the  discreteness  of  data  in  the  spatial 
domain.  If  this  is  not  done,  spatial  aliasing  will  oeear  and  will  be 
•ppnroat  in  the  oatpat.  The  regaireaeat  which  we  iapose  is  that  the 
transverse  component  of  the  wave  vector  in  the  phase  of  the  x  integral  asst 


be  bounded  by  the  tranverse  spatial  Nyquist  frequency  assoeiated  with  the 
saapling  rate  or  spacing  of  the  traces,  as  discussed  next. 

Let  os  denote  the  trace  separation  in  the  transverse  by  Ax.  The 
transverse  wave  vector  can  be  determined  Iron  the  phase  2mx  by  calculating 
the  partial  derivative  of  that  function  with  respect  to  x.  This  derivative 
is  equal  to  the  quotient  of  derivatives  of  2ot  and  |x  -  (|  —  with  t  and 
|x  -  l |  defined  in  (2S)  —  with  respect  to  p  at  fixed  f.  That  ia. 


2w 


dr 

Tv 


-  2m 


P 

cTUT 


(26) 


Ve  introduce  the  frequency  in  Us.  f  ■  «/2x,  and  require  that  the  transverse 
wave  nuaiber  satisfy  our  Nyquist  bound  even  for  the  naxinun  frequency*  any 

f  .  Then 

♦ 


2*2*  f  p  2 

sm — i-n:  • 


(27) 


the  condition 
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p 


TZ  * 


(28) 


This  i»  tli*  additional  constraint  on  p  which,  in  tarn,  impose ■  an  additional 
conatraint  on  the  limits  of  integration  (summation)  in  z. 

The  next  qnestion  which  arises  is  how  to  estimate  the  "jump*  in  a  at  an 
interface  in  terms  of  the  reflectivity,  8.  This  is  discussed  in  Cohen  and 
Bleistein  (1979b]  and  Bleistein  [1984]  for  the  case  of  a  •box*  filter.  Let 
as  define  Aa  as  the  increment  in  a  at  snoh  an  interface  and  P  as  the  peak 
vain*  of  8  *1  an  interface  location  on  an  output  trace.  The  theory  then 
predicts  that 


An  *  P 


c(z) 

4<f+ -  I_>  * 


(29) 


In  this  equation,  f_  is  the  minimum  frequency. 

For  a  more  general  filter,  we  must  replace  the  difference  of 
frequencies  by  the  area  under  the  filter  functional  that  is,  if  A  denotes 
the  area  under  the  filter,  then 


Aa  •  P 


c(s) 


(SO) 
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Finally,  there  ie  tone  poet -process ins  possible  to  oonpensete  partially 
for  the  linearity  of  the  uderlyiag  theory.  This  work  was  reported  in  Begin 
and  Cohen  [1984].  Ve  quote  here  the  reanlts  of  one  iaportant  post-process 
and  refer  the  reader  to  that  paper  for  farther  detail. 

In  order  to  obtain  an  estiaate  of  the  velocity  increaeat  across  an 
interface,  one  uses  the  linearised  estiaate  dednced  froa  (1), 


Ac  .  ■  -1/2  e(s)  An  . 

eat 


(31) 


Here,  Ac  denotes  the  estiaate  obtained  aaaerieslly.  Begin  and  Cohen  then 
eat 

snggeat  as  an  iaproved  estiaate  of  the  velocity  increaeat 


2e(s)  Ac 
_ _ cat 

2c(s)  -  Ac 

cat 


(32) 


Ve  reaark  that  all  each  post-processing  refiaeaeats  arc  only  worthwhile 
if  the  original  data  is  anch  that  there  is  reason  to  believe  that  at  least 
relative  tree  aaplitade  was  preserved.  Vhea  this  is  not  the  ease,  each 
post-processing  is  aaaeeessary.  Iapleacatation  of  (24)  sad  (25)  still  aey 
be  eaployed  to  yield  a  strootnral  Inversion  or  aigratioa  of  the  tiae  logs. 
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«•  shall  bov  discuss  CPU  tin  for  the  iapltatmtatioa  of  the  five  stops 


f. 


listed  ia  the  previous  section.  Let  as  define  the  following: 


The  number  of  points  ia  the  transverse  et  which  the  output  is  to  be 
eelealsted. 


N^.  The  nnsiber  of  points  in  depth  et  which  the  output  ia  to  be  eelculeted. 

The  somber  of  tine  points  on  each  input  trace. 

N  The  amber  of  traces, 
s 

It  should  be  noted  here  that,  ia  praetioe.  N^.  is  saaller  than  N^.  The 
reason  for  this  is  that  the  s sapling  rate  determines  a  Nyquist  frequency 
which  is  usually  larger  than  the  maximum  frequency  of  usable  data.  Thus, 
there  is  a  limit  to  the  resolution  of  the  output  based  upon  this  maximum 
usable  frequency  and  not  on  the  sampling  rate. 

In  our  research,  we  have  concluded  that  the  density  of  the  output  meed 
be  no  greater  than  to  produce  four  points  on  the  main  lobe  of  the 
baadlimited  delta  functions  which  are  being  depicted.  In  Appendix  P,  this 
analysis  is  discussed,  then  those  results  are  implemented  with  authors 
typical  of  geophysical  exploration  data,  we  conclude  that  a  typical  value 
for  the  ratio  is  .Sd. 


Inpleaenta t ion  of  each  step  of  the  previous  section  is  estimated  ss 


follows : 


Step  (i)  0(N  N  log  N  )  to  compute  the  Fourier  transforms. 

x  x  a  x 

Step  (ii)  0(N  N  )  to  multiply  all  of  the  trsasforas  bp  a  function  of 

XT 

frequency  and  then  as  in  Step  (i)  to  invert  the  Fourier 
transform. 


Step  (iii)  0(N  Nj. )  to  develop  the  table  of  p  values. 


Step  (iv)  As  in  Step  (iii). 


Step  (v)  0(N^Nj.NO  to  calenlate  the  final  integral  for  the  output. 

In  Step  (v).  we  use  N’  to  connote  that,  in  fact,  the  aperture  over  whieh  the 
integral  is  to  be  ealcmlated  is  a  fraction  of  the  total  auaber  of  traeea  - 
on  average,  perhaps  one-half. 


Steps  (iii)  aad  (iv)  are  computations  whieh  need  not  be  performed  la 
the  ease  of  constant  c.  Ve  note  that  the  *0*  estimate  for  each  of  these 
steps  oaa  be  sigaifieaatly  smaller  thaa  or  eomparable  to  the  estimate  for 
Step  (v).  His  is  so  because  might  ho  ss  large  as  aad  may  he  a 
small  fraetiea  of  or.  at  most,  equal  to 


at 


7 


V«  use  thru  murieil  examples  both  to  illustrate  the  preceding  theory 
and  to  deaionstrate  the  differenees  between  it  and  other  migration  methods. 
The  first  example  is  relatively  simple i  it  will  be  nsed  to  provide  a 
comparison  of  imaging  methods.  The  second  and  third  are  more  complicated 
and  realiatic. 

Pig.  2(a)  shows  a  synthetic  time  section  from  a  dipping  planar  structure 
in  a  region  where  the  velocity  varies  linearly  with  depthi  the  dip  angle  is 
60  degrees  and  the  velocity  function  is 

c ( x)  -  (5000  +  4x)  ft/sec.  (33) 

The  trace  spacing  is  50  ft.  Vith  c(x)  taken  as  the  reference  velocity,  the 
reconstruction  as  described  in  Section  5  is  shown  in  Pig.  2(b).  The  time 
section  has  been  eorreetly  migrated*  the  apparent  loss  of  high  frequency 
content  in  the  migration  is  due  to  imposing  the  constraint  (28).  Por 

comparison,  we  show  in  Pig.  2(c)  the  results  of  a  Kirchhoff  migration  of  the 
same  time  section.  Straight  ray  paths  were  assumed  for  this  migration,  and 
the  travel  times  for  diffraction  curves  were  computed  by  using  the  rms 

velocity  (converted  from  a  function  of  time  to  a  function  of  depth),  as 
discussed  by  Schneider  (1978).  Prom  this  extreme  example,  which  features  a 
steeply  dipping  event  and  large  background  velocity  variations,  one  oan 

easily  sec  the  advantages  of  computing  trsveltimcs  along  the  curved  re y 

patha  defined  implicitly  in  Sq.  (25).  The  firehkoff  output  is  incorrectly 
pieced  and  the  planar  reflector  is  depicted  es  s  curved  surface. 
Additionally,  a  k-f  migration  was  performed  on  this  data  after  'stretching* 
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the  tiae  axis  [Stolt,  1979].  This  Migration  was  totally  unsuccessful  in 
iaagiag  the  reflector.  This  ia  becanse.  after  the  tiae  axis  was  stretched, 
the  reflection  event  on  the  resnlting  tiae  section  had  an  apparent  dip 
greater  than  45  degrees,  asking  the  aigration  equation 


sin(aigrated  dip)  ■  tan(apparent  dip) 


(34) 


Meaningless. 

The  second  exaaple  (Pig.  3)  shows  a  salt  doae  in  a  aediua  whose 
layering  is  otherwise  nearly  horixontal.  Because  the  objective  is  to  iaage 
the  flanks  of  the  salt  done  and  not  its  interior,  aigration  with  a  reference 
velocity  which  varies  with  depth  only  is  appropriate.  Fig.  3(a)  shows  the 
aodel  and  Fig.  3(b)  shows  a  finite-difference  tiae  section  froa  the  aodel. 
Part  of  the  salt  doae  is  overhung:  its  dip  exceeds  90  degrees.  The  aodel 
velocity  varies  froa  6000  ft/sec  at  the  top  of  the  aodel  to  16000  ft/sec  at 
the  bottoai  it  ia  piecewise  constant.  The  trace  spacing  is  100  ft.  The 
reconstruction  in  Fig.  3(c)  shows  the  eorreet  iaaging  of  dips  op  to,  but  not 
beyond,  vertioal.  For  a  coaparison  with  a  reverse-tine  finite  difference 
aigration  [Vhitaore  1978],  see  Fig.  3(d)  froa  Vhitaore  [1983].  The  finite 
difference  aigration  eorreotly  iaage s  dips  exceeding  vertioal,  but  its 
iapleaentation  is  considerably  aore  tiae-consuaing  than  the  approach  studied 
here. 
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The  third  exaaple  also  ooaaists  of  a  aalt  doaa  in  a  Bedim  whose 
layering  ia  otherwise  Marly  horiaoatal.  The  aodel  is  shown  in  Fig.  4(a). 
along  with  aeleeted  ray  paths  froa  refleeting  interfaces  to  the  npper 
surface.  The  aodel  Telocity  varies  froa  S2S0  ft/sec  at  the  top  of  the  aodel 
to  14100  ft/see  at  the  bottoa.  and  the  trace  spacing  is  2S  ft.  We  note  that 
the  ray  paths  are  piecewise  straight  line  segaents,  indicating  that  the 
aodel  consiata  of  piecewiae  constant-velocity  layers,  contrary  to  our 
specification  of  a  continuous,  pieoewise  linear,  depth  dependent  velocity 
layering.  Nevertheless,  the  aigration  (Fig.  4(b))  was  successful.  The  aost 
notable  errors  are: 

(1)  The  inaccurate  aapping  of  the  interface  directly  below  the  salt  done. 
This  occurred  because  the  acoustic  velocity  of  the  salt  is  higher  than 
that  of  the  sediaents,  leading  to  reduced  traveltiaes  for  reflection 
events  whose  rays  have  travelled  large  diatanoes  through  the  salt.  The 
one-diaeasioaal  reference  velocity  failed  to  aeeouat  for  this  fact. 

(2)  The  appearance  of  diffraction  sailes  on  the  aigrated  section.  These 
are  due  to  the  artificial  truncation  of  reflected  energy  on  the  tiae 
sectioni  that  ia,  the  tiae  section  does  not  contain  diffraction 
effeeta.  This  is,  therefore,  a  liaitation  of  the  synthetic  tiae  data 
and  not  of  the  aigration.  In  fact,  any  aigration  aethod  will  correctly 
interpret  those  truncated  events  as  nearly  circular  reflector 
coat inuat ions.  Finally,  Fig.  4(d)  shows  the  effects  of  aodifying  the 
aati-spatial-aliasiag  constraint  (28).  With  the  4  in  the  denoainator 
replaced  by  8,  spatial  aliasing  ia  further  liaited  along  with  the 
ability  to  iaage  the  aost  steeply-dipping  events  (whieh  appear  on  Fig. 
4(a)  as  those  aost  lihely  to  be  aliased). 
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Coaputer  run  tiaes  for  these  end  other  experiaents  indicate  that  the 
speed  of  this  algorithm  is  about  the  saae  as  that  of  a  coaparably  coded 
(FORTRAN)  k-f  aigration  routine. 
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V*  km  ntndid  tk<  isltidlataiioul  Bora  imnioa  foraalim  tad 
algorithai  to  tk*  oaae  vkm  tki  kaowa  roforoaoo  velocity.  fxoa  which 
portarbatioaa  m  to  ho  competed,  art  a  faaotioa  of  depth.  The  aaalytieal 
rooalt  has  aada  mac  of  a  liaoariiatioa  (aaall  portarbatioaa  froa  the 
roforoaoo  velocity)  aad  a  high  froqaeaey  assaaptioa  (that  aoiaaic  data  aro 
froqaoaey  baadliaited,  with  area  the  lowoat  awailablo  froqmoacioa  coaaidorod 
to  bo  "high"  for  the  probloa  at  haad).  Tho  rodaotioa  to  aa  iaploaoatablo 
algoritha  has  aado  farther  mao  of  tho  high  froqaoaey  aaeaaptioa,  roamltiag 
ia  aa  offioioat  atraetaral  iaworaioa  aothod. 


Vo  wioh  to  thaak  Daa  Vhitaore  for  prowidiag  tho  fiaito-difforoaoo 
ayathotie  data*  aad  laoeo  Prodaetioa  Coapaay  for  poraiaaioa  to  pmbliah  theae 
rooaitt. 

Thia  xoaoaroh  projoet  wao  iaitiatod  while  tho  first  aathor  was  a 
wi alter  at  tho  ioooo  leaearoh  Coator  la  Talsa  dariag  tho  aoatho  of  Febraary 
aad  Varoh*  IMS.  This  aathor  wlahoa  to  osproaa  hia  deep  appreoiatiatioa  to 
haoe©  for  prowidiag  a  frioadly  aad  stiaalatlag  oawlroaaoat  ia  whioh  to  work. 
There  ia  ao  qmoatioa  ia  ay  aiad  that  tho  qaality  of  that  eawiroaaeat 
ooatribatod  sigaifioaatly  to  tho  smeeosa  of  thia  projoet. 
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APPENDIX  A:  AMMPIOT1C  DOVNBOINB  SQUrnCM 


The  purpose  of  this  appendix  is  to  derive  the  asymptotic  downgoing 
solution,  (8),  to  equation  (5).  Ve  begin  by  introducing  the  spatial  Fourier 
transfora  of  the  solution,  defined  by 


f 


u(x,z,t») 


expl-iyx}  dx  dy,  ji  = 


(A-l) 


By  applying  this  Fourier  transfora  to  equation  (5),  we  obtain  the  following 
ordinary  differential  equation  for  the  function  X: 


d*S 


dx 


®  +  p*  8  «  -6(x)  exp(-i|i*$)  , 


(A-2) 


p,  (x.jiia) 


c*(x) 


s 

“  M  » 


p*  ♦  p*  < 

rs 


c*(x) 


le  have  only  defined  ^  for  the  ranges  of  its  variables  which  wake  it  real# 
and  we  shall  only  concern  outselves  with  the  asymptotic  solution  in  that 
range  aa  well. 

Asyaptotie  theory  for  ordinary  differential  equations  —  see  for 
esaaple  Codding ton  and  Levinson  [1955]  —  indieates  that  two  linearly 

independent  solutions  to  the  homogeneous  form  of  (A2)  are  given 
asymptotically  by 


M 


a  +  *i  E'iJ 


3+  “ 


T 


-.  #  *  »*» 


M/z'.fcJw)  d* ’ 


(A-3) 


In  this  equation,  a+  are  eonatants.  The  eolation  u+  is  outgoing  for  s  ->• 
and  the  solution  9  is  outgoing  for  z  — Since  (5)  is  homogeneous  for  z 
nonzero,  ve  conclude  that 


8  -  9+,  +z  >  0 


(A-4) 


with  only  the  constanta  a+  to  he  determined. 

In  order  to  detenaine  thse  constants,  we  must  impose  two  conditiona  at 
z  *  0,  those  conditions  characterising  the  distributional  nature  of  the 
source.  The  conditions  are  that  the  function  itaelf  must  be  continuoua  and 
the  first  derivative  of  the  function  must  have  a  •jump*  equal  to  the 
integral  of  the  right  side  on  an  interval  containing  the  origin  in  s.  These 
two  equations  are 


a 


♦ 


a 


♦ 


-  a  _  •  0  , 

enp C-ig* J) 
♦  a_  ■  -  ; 


(A-5) 


after  a  division  of  comiM  factors.  The  aw  of  those  equations  yields  a+ 
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which  is  ths  osl  y  coast  sat  of  iatsrest  ia  oar  farther  disoassioa. 
lowers ioa  thca  yields  eqastioa  (1). 


Poarisr 


IH  B: 


isnmne  nmuav  or  as) 


Tki  fsrpoit  of  tkii  ifpndix  it  to  doriwo  ti*  rot  alt  (14),  wkiek  it 
tko  atyaptotie  exposa ioa  of  tko  ji”i*tegral  im  (IS).  Wo  tkall  bogia  bo 
reminding  tbo  rotdor  of  tko  bttio  aslti-diaoaaioaal  stationary  pkato 
foxaalt.  For  a  derivation,  too  Bloittoia  (1984),  for  oztaplo. 

Lot  at  tappoto  tkat  1(1)  it  t  salti-fold  iatogrtl  of  tko  fora 


1(1) 


f(k)  ozp(ili(k))  dki*”dk> 
D 


(B-l) 


la  tkit  oqaatioa,  D  doaotoi  aoao  doaaia  la  k  -  tpaeo.  Tko  iatogral  (IS)  it 
of  tkit  font  ia  tko  wariabloa  k  witk  a  -  2.  la  a  aSbaoqaoat  appendix,  wo 
aka  11  kawo  aood  of  tkit  roaal t  witk  a  •  4.  Oar  iatoatloa  it  to  atato  tko 
formal  a  for  tko  loadiag  ton  of  tko  aeymptotio  oapaatioa  of  tkit  iatogral 
for  "large"  waleet  of  tko  paroaotor,  1.  At  ozplaiaod  ia  tko  tost,  wo  kawo 
aot  root  at  tko  iatogral  ia  a  fon  la  wklok  1  it  oxplioit.  Tkoroforo,  wo 
tkall  apply  tko  foraala  bolow  witk  1  ■  1. 

Lot  at  tappoto  tkat,  la  D,  k  la  tko  oaly  tlaplo  ttatioaary  poiat, 

”o 

tkat  it, 

W*(*o>  -  0  .  (B-l) 


bat  tko  aatria  witk  oloaoata 


-  If 


i.j  -l.***.n 


(B-3) 


3*f(k  ) 

fij  "  ak^  * 


is  nonsingnlar .  Then.  the  leading  ten  of  tks  asymptotic  expansion  of  tka 
integral  (B-l)  is 


I,  la/ 2  f(k  ) 

t— -  I  - — — -  exp(iht(h  )  +  ix(sga  iA1)/4]  (B-4) 

1  J  11 


la  tkis  eqaation,  sga  <  is  the  signs tare  of  the  aatrix  f  .  The  sigastare 

*  J 

is  the  nnaber  of  positive  eigeavalaes  aiaas  the  nnaber  of  negative 
eigeavalaes  of  the  aatrix.  If  the  doaain  0  has  aore  stationary  points,  one 
siaply  adds  «p  the  eoatribntions  froa  eaeh  of  thea. 


Let  as  now  consider  the  integral  (13)  and  set 


f(g)  *  sgn  «[#(*, gi«)  ♦  #(s,2k  ♦  »  (B-5) 

with  f  defined  by  (9).  Differentiation  with  respect  to  g  yields  the 
following  reaalts: 


40  - 


(B**7) 


In  these  •qutiou,  both  indices  take  oa  thi  valnes  1  and  2  and  6^  denote* 
tha  Croat eker  dalta  fvaetioa. 

Ha  phase  has  a  stationary  point  whan  i  -  1,2.  Tha  wain* 
of  tha  phasa  at  tha  stationary  point  provides  tha  phase  in  tha  integral  in 
(14) .  Farthemore , 


and  arrive  at  (14). 


Ia  this  appendix,  we  shall  derive  the  result  (19)  froa  the  integral 
(17).  To  do  this,  we  shall  first  calculate  the  leading  ten  of  the 
asymptotic  expansion  of  the  fourfold  integral  in  the  variables  k  and  x.  In 
order  to  iapleaent  the  stationary  phase  fonula  (B-4).  we  think  of  the 
latter  two  variables  as  being  k  and  k  . 

S  4 

Let  us  define  the  phase  to  be  analysed  to  be 

<  »  2(sgn  a)  I  |i  (x'.kjw)  dx'  +  2k*(£  -  x)  .  (C-l) 

Jr 


Ve  take  the  derivative  of  this  phase  function: 


(C-2) 


The  eoaditioaa  that  the  phase  be  stationary  yield  the  solution 


xi *  «i 


i  -  1.2 


(C-3) 


0. 


Ve  aust  now  calculate  the  Matrix  of  aeoond  derivative*  at  the 
atationary  point.  This  is  greatly  facilitated  by  eaploying  (C-3)  in  the 
calculation  process  to  disregard  contributions  which  are  aero  at  the 
stationary  point.  Ve  find  that  at  the  stationary  point. 


< 


11 


"is 


c(z')  da* 

r 


< 


it 


(C-4) 


with  all  other  tens  being  zero.  For  this  Matrix,  we  find  that 


dat  -  16,  sgn  lij  -  0 


(C-5) 


and 


f 


■  2a 


da* 

cTTT  * 

r 


(C-6) 
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4* 


I 


L  J 


o(x')  dx* 


The  w  integration  can  now  be  carried  oat  to  yield 


c(f)  c(0) 


a({. z)  sjjj-  dz'/c(z’)j 

c(x)  I  c(z')  dx' 

J0 


Ve  now  apply  to  (C8)  the  rale. 


I6(f(x))  g(z)  dx  »  -fir? —  , 

|r<M 


(C-8 


(C-9 


with  f  being  the  only  (staple)  xero  of  f(x)  in  the  doaain  of  integration 
This  yields  the  resnlt  (19). 


Ammix  D:  0QHPAKI8QN  WIT*  BAILIES  KESULT 


In  this  appendix  we  shall  disease  the  coaparison  of  (23)  with  the 
corresponding  resalt  froa  Cohen  end  Bleistein  [1979s].  In  the  letter,  the 
foraale  for  a,  in  the  two  end  one-helf  diaensionel  case,  is  stated  in 
equations  (9)  and  (10).  Ve  quote  that  result  with  changes  in  notation  to 
correspond  to  the  notation  here:  we  aust  interchange  the  roles  of  x.  (  and 


x.  The  foraula  for  a  then  becoaes 


o(t.t)  - 


CD  f 


(D-l) 


•  r(r  -  t)  Og(x.O.t)  exp(2i[ki((  -  x)  -  k^f]  ♦  iat]  . 


In  this  equation. 


-  o.gnk  lkx  ♦  k. 


(D-2) 


Ve  ahull  now  siaplify  the  t iae  doaaia  integrals.  To  begin,  let  ua 


-  43  - 


define 


I  ■  f  dx  f  dtx(x  -  t)  D_(x,0,t)  exp(iwx)  .  (D-3) 
J0  J0  8 

la  tbit  integral,  we  with  to  interchange  the  order  of  integration.  In  doing 
to,  we  take  the  point  of  view  that  m  hat  a  positive  imaginary  part.  Indeed, 
there  it  jnttif ication  for  thit.  The  original  problen  waa  "causalj*  that 
it,  the  experinentt  were  carried  ont  etarting  from  tone  finite  time.  That, 
the  Foarier  trantform  to  frequency  domain  it  originally  defined  and  analytic 
in  tome  upper  half  w-plane.  Before  analytically  continuing  the  trantform 
down  to  the  real  axit,  the  interchange  it  justified.  After  carrying  out  the 
interchange,  we  obtained 


I 


I  dt  U_(x,0,t)  I  x  dx(x  -  t)  exp(iwx)  . 
0S  Jt 


(D-4) 


Ve  now  calculate  the  x  integral  by  integrating  by  parte  and  keeping 
the  leading  order  term  in  «  for  large  |«e|.  The  result  it 


I 


-  I  dt  t  Ug(x,0,t)  exp(iwt)  • 
*  ^0 


(D-5) 
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which  we  ei 


1 


i 


1  d 


(D-6) 


I 


.  *  dw 
iu> 


I  dt  0s<x,0,t)  exp(iwt} 
^0 


We  substitute  this  result  into  (D-l)  and  also  introduce  the  change  of 
variable  of  integration  froa  k^  to  w  via  the  equation 


k 

t 


sgn  w 


(D-7) 


and  obtain 

0(5, f)  *  -jp  J  dx  j*  dk^  exp{2i[k1({  -  x)  -  f  1 } 

Oft 

•  -jj- J  Ug(x,0.t)  expUwtJdt  . 


(D-8) 


The  theory  developed  in  Cohen  and  Bleisiein  [1979a]  indicates  that  «e 
obtain  the  reflectivity  function  by  aultiplying  the  integrand  by  2ii*/c, 
which  yielda  the  following  result: 


-  ~  J  4*  J  dk4  |  dto  exp(2i[k1((  -  x)  -  k.tl) 

*  |  0g(x,0,t)  exp(itot)  dt 

An  integration  by  parts  in  m  now  yields  the  result 


(D-9) 
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exp{2i[kx((  -  x)  -  kJf]} 


(D-10) 

00 

*  I  dt  Dg(x,0,t)  exp(iwt}dt 
■^0 


When  (23)  is  specialized  to  the  case  of  constant  e,  the  result  is  exactly 
the  sane  as  (D-10)  when  one  Bakes  the  identification 


|ig  sgn  *  -  k,  . 


(D-ll) 


tfmn  B:  Asnrronc  expansion 


We  shell  describe  here  the  asynptotic  expansion  bp  the  Method  of 
stationery  phase  of  the  integral  with  respect  to  k  in  (23).  We  begi?  by 
rewriting  that  inegral  as 


nc*( f)c(O) 


»  i 

| is ) dee  dx  dt  Ug(x.O,t) 

4  J  * 


» 

dk  exp(2i<  +  iwt) 

|i~  (s.kiM) 

i 


(B-l) 


In  this  equation, 

f 

f  ■  k((  -  x)  -  agnw  I  |i|(z.kt«)  dx.  (B-2) 

To  calcnlate  this  integral  by  the  Method  of  stationary  phase,  we  need  the 
first  and  second  derivatiwes  of  the  phase  function.  They  are  given  by 


Sf 

Ik 


{  -  i  +  k  agn  « 


dx 

7T7E757 


(B-3) 


and 
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d*| 


1 


dz 


i .  .  t 


(E-4) 


dk  I  c  (z)(it  (z.kiw) 


The  condition  that  the  phase  be  stationary,  that  is,  that  df/dk  *  0,  is 


x  -  5  =  k  sgn 


"  A 


(E-5) 


Ve  set 


k 


7-  p  sgn  (x  -  {) 


(E-6) 


and  then  (E-5)  is  just  the  ray  equation,  which  is  the  first  equation  in 
(25).  Furthermore,  at  the  stationary  point, 

|i,(z.k,M)  -  J  n*(s)  -  p*  ,  I  »  irc(p,f)  , 

(E-7) 

-S 

— -  ■  c(0)  m*  S(f,p)  sgn  u  . 
dk 

Substitution  of  these  values  and  application  of  the  stationary  phase  formula 
to  (E-i)  yields  the  result.  (24),  (25). 
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tflBBB  F:  SUPLIM  BATS  OP  OUTPUT  IN  BPPTB 

We  »h»ll  diiesit  here  the  question  of  the  stapling  rate  in  depth  of  the 
ontpnt  of  i  seisaic  inversion  algorithm.  These  considerations  are  not 
peenliar  to  inversion)  they  apply  to  any  siailar  imaging  technique. 

Let  us  consider  the  final  step  of  the  integral  process  to  iaage  a 
reflector.  For  siaplicity  ve  deal  with  the  case  of  a  horizontal  reflector 
in  one  spatial  diaension.  Then,  the  last  integral  to  be  perforaed  is  of  the 
fora 


I(z)  ■  I  ezpC4nifz/c)  df 
"B 


(F-l) 


In  this  equation,  the  band  B  over  which  the  integral  is  to  be  calculated  is 
a  syaaetric  pair  of  intervals,  <-f+,  -f_),  (f_,  f+).  Calculating  this 
integral  explicitly  and  using  a  trigonoaetric  identity  for  the  difference  of 
sines  of  two  angles  yields  the  result 


I<*> 


c  cos [2n(f+  +  f_)z/cj  sin[2x(f+  -  f_)z/oJ 
' .  xz 


(F-2) 


The  zeroes  of  this  function  nearest  to  the  origin  oeeur  when  the  eosine 
factor  in  this  equation  is  zero  or  when  the  arguaent  of  the  eosine  factor  is 
x/2.  The  distance  between  these  two  zeroes  is  deteraincd  by  setting  the 
erguaent  of  the  cosine  equal  to  n. 
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Oa  empirical  ground*,  we  require  that  the  sampling  rate  ia  depth  be 
aaoh  that  there  are  four  sample  poiats  ia  this  interval.  Therefore,  if  Ax 
is  the  sample  interval, 

2a (f  +  f  )  4  Ax 

- - -  „  »  (F— 3 ) 

o 

or. 


Ax  “ 


8(f+  +  f_ 


(P-4) 


Let  us  defiae  the  maximum  depth  to  be  Z  aad  the  BMximum  time 
corresponding  to  that  depth,  to  be  T.  Then,  aeglectiag  the  practical 
aspects  of  processiag  to  the  maximum  depth,  we  set 


Nx  "  IT  "  4At(f+  +  f->  Nt 


(P-5) 


Por  a  sample  rate  of  4  mils  aad  a  useful  frequency  range  of  5-30  Us,  we 

obtain  from  (P6)  the  estimate.  N  /N  -  .56. 

x  t 


Figure  1.  The  bandliaited  singular  function  of  n  surface. 

Figure  2(a).  Synthetic  tiae  section:  dipping  planar  refleotor  in  a  aediua 
where  the  velocity  increases  linearly  with  depth.  Dip  angle  is  60  degreesi 
velocity  is  (S000  +  4s)  ft/sec. 

Figure  2(h).  Depth  section  reconstructed  via  the  Born  procedure.  The 
reflector  has  been  properly  located. 

Figure  2(c).  Depth  section  reconstructed  via  Kirehhoff  aigration  using  ras 
velocities.  Because  of  errors  in  coaputing  traveltiaes  for  diffraction 
curves,  the  result  is  a  severe  overaigration. 

Figure  3(a).  Salt  doae  aodel. 

Figure  3(b).  Synthetic  tiae  section. 

Figure  3(c).  Deconstructed  depth  section.  All  dips  up  to  90  degrees  have 
been  properly  aigrated. 

Figure  3(d)  Output  of  WMtaore’s  reverse-tine  finite  difference  aigration. 
(With  the  peraission  of  the  author.) 

Figure  4(a).  Salt  doae  aodel  with  ray  paths  froa  one  reflecting  horizon. 
Figure  4(b).  Synthetic  tiae  section. 

Figure  4(c).  Keconstrueted  depth  section. 

Figure  4(d).  Keconstrueted  depth  section*  the  anti-aliasing  constraint  (28) 
has  been  modified. 
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AMBICr 

Ik*  purpose  of  this  pspsr  is  to  doseribo  an  extension  of  the 
anltidisensionl  Born  inversion  technique  [Cohen  end  Bleistein,  1979s]  for 
eooestie  wares.  In  that  earlier  work,  a  perturbation  in  reference 
soundspeed  was  deterained  bp  assuming  that  the  reference  or  background  speed 
was  constant.  In  this  extension,  we  allow  the  reference  speed  to  be  a 
function  of  the  depth  variable,  s.  bat  still  require  that  it  be  independent 
of  the  transverse  variables.  The  output  of  this  sethod  is  a  high  frequency 
bandlinited  reflectivity  function  of  the  subsurface.  The  reflectivity 
function  is  an  array  of  bandlinited  singular  functions  scaled  by  the  normal 
reflection  strength.  Bach  singular  function  is  a  Dirac  delta  function  of  a 
scalar  argaswnt  which  measures  distance  normal  to  a  reflecting  interface. 
Thus,  the  reflectivity  function  is  an  indicator  nap  of  subsurface  reflectors 
which  is  equivalent  to  the  map  produced  by  migration.  In  addition  to  the 
assumption  of  small  perturbation,  the  method  requries  the  assumption  that 
the  reflection  data  reside  in  the  high  frequency  regime,  in  a  well-defined 
sense . 

The  method  is  based  on  the  derivation  of  an  integral  equation  for  the 
perturbation  in  soundspeed  from  a  known  reference  speed.  When  the  reference 
speed  is  constant,  the  integral  equation  admits  an  analytic  solution  as  a 
multifold  integral  of  the  reflection  data..  Further,  high  frequency 
asymptotic  analysis  simplifies  this  integral  considerably  and  leads  to  sa 
extremely  efficient  numerical  algorithm  for  computing  the  relfectivity 
function.  In  a  paper  by  Bleistein,  Cohen  and  lagia  [1984],  the  development 
of  a  computer  code  to  implement  this  constant  reference  speed  solution  is 
described. 

For  a  depth  dependent  reference  speed,  e(s),  we  can  no  longer  invert 
the  integral  equation  exactly,  lowever,  we  earn  write  down  an  asymptotic 
high  frequency  approximation  for  the  kernel  of  the  integral  equation  and  aa 
asymptotic  eolation  for  the  perturbation.  The  computer  implementation  of 
this  result  is  designed  along  the  same  lines  as  the  code  for  the  constant 
background  ease.  In  tests  of  processing  time,  we  find  that,  at  worst,  the 
total  processing  time  for  this  algorithm  with  depth  dependent  background 
soundspeed  is  about  the  same  as  for  a  comparably  programmed  k-f  algoritha 
with  constant  background .  By  worst  we  mean  that  we  ehoose  the  aspect  ratic 
—  the  number  of  traees  divided  by  the  numbs*  of  points  per  trues  —  to  be 
optimal  for  the  k-f  algorithm.  Vs  present  examples  which  desMustrste  the 
method  isqtlementsd  as  a  migration  technique  and  compare  with  the  applieatioa 
of  alternative  migration  algorithms.  The  examples  we  sheas  were  ones  is 
which  the  objective  is  to  image  the  flanks  of  a  salt  demo. 
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